Summary
This is a preparatory phase to convert the satellite data of elevation details into the format that can be processed by programs and tools in standard languages.
STEP1: Convert GeoTiff files to NumPy files
#mkdir /opt/appfs/bins
python3 /opt/gistools/convert_geotif_to_bin.py --input_path /opt/appfs/ASTER_DEM_TIF_Files --output_path /opt/appfs/ASTER_DEM_Bin_Files --file_regex "AST.*_dem.tif"
<<This will generate ".bin" and ".json" files in output_path, /opt/appfs/ASTER_DEM_Bin_Files. For e-g., ASTGTMV003_N31E068_dem.tif.bin, and ASTGTMV003_N31E068_dem.tif.json>>
For part of tiles (e-g., N28,N29):
python3 /opt/gistools/convert_geotif_to_bin.py --input_path /opt/appfs --output_path /opt/appfs/bins --file_regex "ASTGTMV003_N2[89].*_dem.tif"
Example reviewing code:
import rasterio
python3
fp = r'ASTGTMV003_N35E097_dem.tif'
dataset = rasterio.open(fp)
# Get all metadata
dataset.profile
#<<{'driver': 'GTiff', 'dtype': 'int16', 'nodata': None, 'width': 3601, 'height': 3601, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.000277777777777778, 0.0, 96.9998611111111,
# 0.0, -0.000277777777777778, 36.0001388888889), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'compress': 'lzw', 'interleave': 'band'}>>
dataset.meta
#<<{'driver': 'GTiff', 'dtype': 'int16', 'nodata': None, 'width': 3601, 'height': 3601, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.000277777777777778, 0.0, 96.9998611111111,
# 0.0, -0.000277777777777778, 36.0001388888889)}>>
# No. of Bands in tif
print(dataset.count)
#<<1>>
# Image resolution
print(dataset.height, dataset.width)
#<<3601,3601>>
# Coordinate Reference System
print(dataset.crs)
#<<EPSG:4326>>
# Get lat/long boundingbox
print(dataset.bounds)
#<<BoundingBox(left=96.9998611111111, bottom=34.99986111111112, right=98.00013888888888, top=36.0001388888889)>>
# Get spatial coordinates of pixel at the center of image
dataset.xy(dataset.height // 2, dataset.width // 2)
#<<(97.49999999999999, 35.50000000000001)>>
# Get daratypes of image data
dataset.dtypes
#<<('int16',)>>
# Access image data of first band
band1 = dataset.read(1)
#<<Returns numpy ndarray>>
#array([[4369, 4380, 4388, ..., 3428, 3427, 3430],
# [4392, 4403, 4413, ..., 3431, 3432, 3434],
# [4406, 4413, 4420, ..., 3428, 3429, 3429],
# ...,
# [4371, 4366, 4361, ..., 4237, 4237, 4236],
# [4370, 4362, 4360, ..., 4236, 4236, 4237],
# [4366, 4360, 4355, ..., 4234, 4236, 4238]], dtype=int16)
band1.shape
#<<(3601, 3601)>>
band1[0,0]
#<<4369>>
band1[3600, 3600]
#<<4238>>
# Get resolution
dataset.res
(0.000277777777777778, 0.000277777777777778)
#move one unit right and one unit down in geo space
x, y = (dataset.bounds.left + 0.000277777777777778, dataset.bounds.top - 0.000277777777777778)
# get spatial coord correspondence in band array
row, col = dataset.index(x, y)
#<<(1, 1)>>
band1[row, col]
<<4403>> import rasterio
fp = r'ASTGTMV003_N35E097_dem.tif'
dataset = rasterio.open(fp)
# Get all metadata
dataset.profile
#<<{'driver': 'GTiff', 'dtype': 'int16', 'nodata': None, 'width': 3601, 'height': 3601, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.000277777777777778, 0.0, 96.9998611111111,
# 0.0, -0.000277777777777778, 36.0001388888889), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'compress': 'lzw', 'interleave': 'band'}>>
dataset.meta
#<<{'driver': 'GTiff', 'dtype': 'int16', 'nodata': None, 'width': 3601, 'height': 3601, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.000277777777777778, 0.0, 96.9998611111111,
# 0.0, -0.000277777777777778, 36.0001388888889)}>>
# No. of Bands in tif
print(dataset.count)
#<<1>>
# Image resolution
print(dataset.height, dataset.width)
#<<3601,3601>>
# Coordinate Reference System
print(dataset.crs)
#<<EPSG:4326>>
# Get lat/long boundingbox
print(dataset.bounds)
#<<BoundingBox(left=96.9998611111111, bottom=34.99986111111112, right=98.00013888888888, top=36.0001388888889)>>
# Get spatial coordinates of pixel at the center of image
dataset.xy(dataset.height // 2, dataset.width // 2)
#<<(97.49999999999999, 35.50000000000001)>>
# Get daratypes of image data
dataset.dtypes
#<<('int16',)>>
# Access image data of first band
band1 = dataset.read(1)
#<<Returns numpy ndarray>>
#array([[4369, 4380, 4388, ..., 3428, 3427, 3430],
# [4392, 4403, 4413, ..., 3431, 3432, 3434],
# [4406, 4413, 4420, ..., 3428, 3429, 3429],
# ...,
# [4371, 4366, 4361, ..., 4237, 4237, 4236],
# [4370, 4362, 4360, ..., 4236, 4236, 4237],
# [4366, 4360, 4355, ..., 4234, 4236, 4238]], dtype=int16)
band1.shape
#<<(3601, 3601)>>
band1[0,0]
#<<4369>>
band1[3600, 3600]
#<<4238>>
# Get resolution
dataset.res
(0.000277777777777778, 0.000277777777777778)
#move one unit right and one unit down in geo space
x, y = (dataset.bounds.left + 0.000277777777777778, dataset.bounds.top - 0.000277777777777778)
# get spatial coord correspondence in band array
row, col = dataset.index(x, y)
#<<(1, 1)>>
band1[row, col]
<<4403>>STEP2: Convert Python Numpy binary files to Java short binary files
Running conversion with 16 threads:
export G2K_DATA_HOME=/opt/appfs
export NUMPYBINCONVERTER_THREADS=16
nohup ${GEOTIFF_APP_HOME}/run_NumpyBinConverterService.sh convertAllNumpyBinFiles ${G2K_DATA_HOME}/ASTER_DEM_TIF_Files ASTGTMV003.* ${G2K_DATA_HOME}/ASTER_DEM_Bin_Files ${NUMPYBINCONVERTER_THREADS} > /opt/geotiff_reader_logs/nohup.log 2>&1 &
<<For 108 tiles, it takes ~10min on 128GB server>>